library(dplyr)
library(INLA)
library(ggplot2)
library(patchwork)
library(inlabru)
# load some libraries to generate nice plots
library(scico)Practical 1
Aim of this practical:
In this first practical we are going to look at some simple models
- A Gaussian model with simulated data
- A Linear mixed model
- A GLM model with random effects
we are going to learn:
- How to fit some commonly used models with
inlabru - How to explore the results
- How to get predictions for missing data points
Linear Model
In this practical we will:
- Simulate Gaussian data
- Learn how to fit a linear model with
inlabru - Generate predictions from the model
Start by loading useful libraries:
As our first example we consider a simple linear regression model with Gaussian observations \[ y_i\sim\mathcal{N}(\mu_i, \sigma^2), \qquad i = 1,\dots,N \]
where \(\sigma^2\) is the observation error, and the mean parameter \(\mu_i\) is linked to the linear predictor (\(\eta_i\)) through an identity function: \[ \eta_i = \mu_i = \beta_0 + \beta_1 x_i \] where \(x_i\) is a covariate and \(\beta_0, \beta_1\) are parameters to be estimated. We assign \(\beta_0\) and \(\beta_1\) a vague Gaussian prior.
To finalize the Bayesian model we assign a \(\text{Gamma}(a,b)\) prior to the precision parameter \(\tau = 1/\sigma^2\) and two independent Gaussian priors with mean \(0\) and precision \(\tau_{\beta}\) to the regression parameters \(\beta_0\) and \(\beta_1\) (we will use the default prior settings in INLA for now).
We can write the linear predictor vector \(\pmb{\eta} = (\eta_1,\dots,\eta_N)\) as
\[ \pmb{\eta} = \pmb{A}\pmb{u} = \pmb{A}_1\pmb{u}_1 + \pmb{A}_2\pmb{u}_2 = \begin{bmatrix} 1 \\ 1\\ \vdots\\ 1 \end{bmatrix} \beta_0 + \begin{bmatrix} x_1 \\ x_2\\ \vdots\\ x_N \end{bmatrix} \beta_1 \]
Our linear predictor consists then of two components: an intercept and a slope.
Simulate example data
First, we simulate data from the model
\[ y_i\sim\mathcal{N}(\eta_i,0.1^2), \ i = 1,\dots,100 \]
with
\[ \eta_i = \beta_0 + \beta_1 x_i \]
where \(\beta_0 = 2\), \(\beta_1 = 0.5\) and the values of the covariate \(x\) are generated from an Uniform(0,1) distribution. The simulated response and covariate data are then saved in a data.frame object.
Simulate Data from a LM
beta = c(2,0.5)
sd_error = 0.1
n = 100
x = rnorm(n)
y = beta[1] + beta[2] * x + rnorm(n, sd = sd_error)
df = data.frame(y = y, x = x) Fitting a linear regression model with inlabru
Defining model components
The model has two parameters to be estimated \(\beta_1\) and \(\beta_2\). We need to define the two corresponding model components.
Note that we have excluded the default Intercept term in the model by typing -1 in the model components. However, inlabru has automatic intercept that can be called by typing Intercept() , which is one of inlabru special names and it is used to define a global intercept, e.g.
cmp = ~ Intercept(1) + beta_1(x, model = "linear")Observation model construction
The next step is to construct the observation model by defining the model likelihood. The most important inputs here are the formula, the family and the data.
The likelihood for the observational model is defined using the bru_obs() function.
Fit the model
We fit the model using the bru() functions which takes as input the components and the observation model:
fit.lm = bru(cmp, lik)Extract results
The summary() function will give access to some basic information about model fit and estimates
summary(fit.lm)
## inlabru version: 2.13.0.9011
## INLA version: 25.09.19
## Components:
## Latent components:
## beta_0: main = linear(1)
## beta_1: main = linear(x)
## Observation models:
## Family: 'gaussian'
## Tag: <No tag>
## Data class: 'data.frame'
## Response class: 'numeric'
## Predictor: y ~ beta_0 + beta_1
## Additive/Linear: TRUE/TRUE
## Used components: effects[beta_0, beta_1], latent[]
## Time used:
## Pre = 1.02, Running = 0.324, Post = 0.0178, Total = 1.36
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## beta_0 1.988 0.009 1.969 1.988 2.006 1.988 0
## beta_1 0.491 0.010 0.472 0.491 0.510 0.491 0
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 118.25 16.72 87.80 117.47
## 0.975quant mode
## Precision for the Gaussian observations 153.24 115.90
##
## Marginal log-Likelihood: 74.36
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')We can see that both the intercept and slope and the error precision are correctly estimated.
Generate model predictions
Now we can take the fitted bru object and use the predict function to produce predictions for \(\mu\) given a new set of values for the model covariates or the original values used for the model fit
new_data = data.frame(x = c(df$x, runif(10)),
y = c(df$y, rep(NA,10)))
pred = predict(fit.lm, new_data, ~ beta_0 + beta_1,
n.samples = 1000)The predict function generate samples from the fitted model. In this case we set the number of samples to 1000.
Code
pred %>% ggplot() +
geom_point(aes(x,y), alpha = 0.3) +
geom_line(aes(x,mean)) +
geom_line(aes(x, q0.025), linetype = "dashed")+
geom_line(aes(x, q0.975), linetype = "dashed")+
xlab("Covariate") + ylab("Observations")Linear Mixed Model
In this practical we will:
- Understand the basic structure of a Linear Mixed Model (LLM)
- Simulate data from a LMM
- Learn how to fit a LMM with
inlabruand predict from the model.
Consider the a simple linear regression model except with the addition that the data that comes in groups. Suppose that we want to include a random effect for each group \(j\) (equivalent to adding a group random intercept). The model is then: \[ y_{ij} = \beta_0 + \beta_1 x_i + u_j + \epsilon_{ij} ~~~ \text{for}~i = 1,\ldots,N~ \text{and}~ j = 1,\ldots,m. \]
Here the random group effect is given by the variable \(u_j \sim \mathcal{N}(0, \tau^{-1}_u)\) with \(\tau_u = 1/\sigma^2_u\) describing the variability between groups (i.e., how much the group means differ from the overall mean). Then, \(\epsilon_j \sim \mathcal{N}(0, \tau^{-1}_\epsilon)\) denotes the residuals of the model and \(\tau_\epsilon = 1/\sigma^2_\epsilon\) captures how much individual observations deviate from their group mean (i.e., variability within group).
The model design matrix for the random effect has one row for each observation (this is equivalent to a random intercept model). The row of the design matrix associated with the \(ij\)-th observation consists of zeros except for the element associated with \(u_j\), which has a one.
\[ \pmb{\eta} = \pmb{A}\pmb{u} = \pmb{A}_1\pmb{u}_1 + \pmb{A}_2\pmb{u}_2 + \pmb{A}_3\pmb{u}_3 \]
Simulate example data
set.seed(12)
beta = c(1.5,1)
sd_error = 1
tau_group = 1
n = 100
n.groups = 5
x = rnorm(n)
v = rnorm(n.groups, sd = tau_group^{-1/2})
y = beta[1] + beta[2] * x + rnorm(n, sd = sd_error) +
rep(v, each = 20)
df = data.frame(y = y, x = x, j = rep(1:5, each = 20)) Note that inlabru expects an integer indexing variable to label the groups.
Code
ggplot(df) +
geom_point(aes(x = x, colour = factor(j), y = y)) +
theme_classic() +
scale_colour_discrete("Group")Fitting a LMM in inlabru
Defining model components and observational model
In order to specify this model we must use the group argument to tell inlabru which variable indexes the groups. The model = "iid" tells INLA that the groups are independent from one another.
# Define model components
cmp = ~ -1 + beta_0(1) + beta_1(x, model = "linear") +
u(j, model = "iid")The group variable is indexed by column j in the dataset. We have chosen to name this component v() to connect with the mathematical notation that we used above.
# Construct likelihood
lik = bru_obs(formula = y ~.,
family = "gaussian",
data = df)Fitting the model
The model can be fitted exactly as in the previous examples by using the bru function with the components and likelihood objects.
fit = bru(cmp, lik)
summary(fit)
## inlabru version: 2.13.0.9011
## INLA version: 25.09.19
## Components:
## Latent components:
## beta_0: main = linear(1)
## beta_1: main = linear(x)
## u: main = iid(j)
## Observation models:
## Family: 'gaussian'
## Tag: <No tag>
## Data class: 'data.frame'
## Response class: 'numeric'
## Predictor: y ~ .
## Additive/Linear: TRUE/TRUE
## Used components: effects[beta_0, beta_1, u], latent[]
## Time used:
## Pre = 0.89, Running = 0.194, Post = 0.0341, Total = 1.12
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## beta_0 2.108 0.438 1.229 2.108 2.986 2.108 0
## beta_1 1.172 0.120 0.936 1.172 1.407 1.172 0
##
## Random effects:
## Name Model
## u IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.995 0.144 0.738 0.986
## Precision for u 1.613 1.060 0.369 1.356
## 0.975quant mode
## Precision for the Gaussian observations 1.30 0.971
## Precision for u 4.35 0.918
##
## Marginal log-Likelihood: -179.93
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Model predictions
To compute model predictions we can create a data.frame containing a range of values of covariate where we want the response to be predicted for each group. Then we simply call the predict function while specifying the model components.
LMM fitted values
# New data
xpred = seq(range(x)[1], range(x)[2], length.out = 100)
j = 1:n.groups
pred_data = expand.grid(x = xpred, j = j)
pred = predict(fit, pred_data, formula = ~ beta_0 + beta_1 + u)
pred %>%
ggplot(aes(x=x,y=mean,color=factor(j)))+
geom_line()+
geom_ribbon(aes(x,ymin = q0.025, ymax= q0.975,fill=factor(j)), alpha = 0.5) +
geom_point(data=df,aes(x=x,y=y,colour=factor(j)))+
facet_wrap(~j)Generalized Linear Model
In this practical we will:
- Simulate non-Gaussian data
- Learn how to fit a generalised linear model with
inlabru - Generate predictions from the model
A generalised linear model allows for the data likelihood to be non-Gaussian. In this example we have a discrete response variable which we model using a Poisson distribution. Thus, we assume that our data \[ y_i \sim \text{Poisson}(\lambda_i) \] with rate parameter \(\lambda_i\) which, using a log link, has associated predictor \[ \eta_i = \log \lambda_i = \beta_0 + \beta_1 x_i \] with parameters \(\beta_0\) and \(\beta_1\), and covariate \(x\).
Simulate example data
This code generates 100 samples of covariate x and data y.
set.seed(123)
n = 100
beta = c(1,1)
x = rnorm(n)
lambda = exp(beta[1] + beta[2] * x)
y = rpois(n, lambda = lambda)
df = data.frame(y = y, x = x) Fitting a GLM in inlabru
Define model components
The predictor here only contains only 2 components (Intercept and Slope).
Define linear predictor
Build observational model
When building the observation model likelihood we must now specify the Poisson likelihood using the family argument (the default link function for this family is the \(\log\) link).
lik = bru_obs(formula = eta,
family = "poisson",
data = df)Fit the model
Once the likelihood object is constructed, fitting the model is exactly the same process, we just need to specify the model components and the observational model, and pass this on to the bru function:
fit_glm = bru(cmp, lik)And model summaries can be viewed using
summary(fit_glm)inlabru version: 2.13.0.9011
INLA version: 25.09.19
Components:
Latent components:
beta_0: main = linear(1)
beta_1: main = linear(x)
Observation models:
Family: 'poisson'
Tag: <No tag>
Data class: 'data.frame'
Response class: 'integer'
Predictor: y ~ beta_0 + beta_1
Additive/Linear: TRUE/TRUE
Used components: effects[beta_0, beta_1], latent[]
Time used:
Pre = 0.825, Running = 0.182, Post = 0.00557, Total = 1.01
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
beta_0 0.915 0.071 0.775 0.915 1.054 0.915 0
beta_1 1.048 0.056 0.938 1.048 1.157 1.048 0
Marginal log-Likelihood: -204.02
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
Generate model predictions
To generate new predictions we must provide a data frame that contains the covariate values for \(x\) at which we want to predict.
This code block generates predictions for the data we used to fit the model (contained in df$x) as well as 10 new covariate values sampled from a uniform distribution runif(10).
# Define new data, set to NA the values for prediction
new_data = data.frame(x = c(df$x, runif(10)),
y = c(df$y, rep(NA,10)))
# Define predictor formula
pred_fml <- ~ exp(beta_0 + beta_1)
# Generate predictions
pred_glm <- predict(fit_glm, new_data, pred_fml)Since we used a log link (which is the default for family = "poisson"), we want to predict the exponential of the predictor. We specify this using a general R expression using the formula syntax.
Note that the predict function will call the component names (i.e. the “labels”) that were decided when defining the model.
Since the component definition is looking for a covariate named \(x\), all we need to provide is a data frame that contains one, and the software does the rest.
pred_glm %>% ggplot() +
geom_point(aes(x,y), alpha = 0.3) +
geom_line(aes(x,mean)) +
geom_ribbon(aes(x = x, ymax = q0.975, ymin = q0.025),fill = "tomato", alpha = 0.3)+
xlab("Covariate") + ylab("Observations (counts)")